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In transport barriers, particularly H-mode edge pedestals, radial scale lengths can 
become comparable to the ion orbit width, causing neoclassical physics to become 
radially nonlocal. In this work, the resulting changes to neoclassical flow and current 
are examined both analytically and numerically. Steep density gradients are con- 
sidered, with scale lengths comparable to the poloidal ion gyroradius, together with 
strong radial electric fields sufficient to electrostatically confine the ions. Attention 
is restricted to relatively weak ion temperature gradients (but permitting arbitrary 
electron temperature gradients), since in this limit a 5f (small departures from a 
Maxwellian distribution) rather than full-/ approach is justified. This assumption is 
in fact consistent with measured inter-ELM H-Mode edge pedestal density and ion 
temperature profiles in many present experiments, and is expected to be increasingly 
valid in future lower coUisionality experiments. In the numerical analysis, the dis- 
tribution function and Rosenbluth potentials are solved for simultaneously, allowing 
use of the exact field term in the linearized Fokker-Planck-Landau collision operator. 
In the pedestal, the parallel and poloidal fiows are found to deviate strongly from 
the best available conventional neoclassical prediction, with large poloidal variation 
of a different form than in the local theory. These predicted effects may be observ- 
able experimentally. In the local limit, the Sauter bootstrap current formulae appear 
accurate at low coUisionality, but they can overestimate the bootstrap current near 
the plateau regime. In the pedestal ordering, ion contributions to the bootstrap and 
Pfirsch-Schliiter currents are also modified. 
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I. INTRODUCTION 



Neoclassical effects in a plasma - the ffows, ffuxes, and currents determined by colli- 
sions in a toroidal equilibrium in the absence of turbulence - set a minimum level of radial 
transportii^. In transport barriers - the pedestal at the edge of an H-mode or internal trans- 
port barriers - neoclassical effects are particularly important for several reasons. First, the 
pressure gradient driven flows and bootstrap current (thought to be determined or at least 
strongly influenced by neoclassical physics even in the presence of turbulence) become large 
due to the small radial scale-lengths. Second, turbulent radial transport is reduced, so neo- 
classical radial transport becomes more relevant. Both the flows and bootstrap current will 
affect the global stability of the transport barrier region. For example, to predict whether 
given plasma proflles are stable to Edge Localized Modes (ELMs) and to predict the nature 
of such ELMs^, accurate calculation of the bootstrap current is essential. 

However, conventional neoclassical calculations are not formally valid in the pedestal. The 
reason is that in conventional neoclassical calculations, the main ion distribution function /i 
is expanded in an asymptotic seriesii^ fi = /Mi + /i + - • • with fi/ fui ^ Pe/r± ^ I, where /mi 
is a Maxwellian, pg = [B / Be)v\/VL is the poloidal ion gyroradius, B = \B\ is the magnitude 
of the magnetic fleld, Bg is the poloidal magnetic fleld, r± ~ |Vlnpi|^^ ~ |VlnTi|~^ is the 
scale- length of ion pressure Pi or ion temperature T-^, vi = \/TLjrni is the ion thermal speed, 
r2 = ZeB j (mic) is the gyrofrequency, Z is the ion charge in units of the proton charge e, mi 
is the ion mass, and c is the speed of light. In the pedestal, is observed to be comparable 
to pq in present experiments. In particular, the density gradient scale length is generally 
comparable to pg. (For this discussion it does not matter whether actually scales with pg.) 
The flrst two terms in the asymptotic series /mi and j\ are then of comparable magnitude, 
so the asymptotic approach breaks down. In the conventional case, the orbit width (~ pe) is 
thin compared to the equilibrium proflles, so neoclassical effects are radially local: the flows 
on a given flux surface depend only on the physical quantities and their radial gradients at 
that surface. However, in a transport barrier where the ion orbit width is not small relative 
to the equilibrium scales, the ions will sample a range of densities and temperatures during 
their orbits. Accordingly, ion flows on a given flux surface are influenced by equilibrium 
parameters from neighboring flux surfaces that lie roughly within a poloidal gyroradius. 
Thus a radially global (i.e. nonlocal) calculation is required for the ion physics. A nonlocal 
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calculation is unnecessary for electrons since their orbit widths are \fm^m\ times smaller 
than ion orbit widths, but the electron distribution is nonetheless modified due to collisions 
with the modified ion distribution^. 

In the conventional local theory, a natural scale separation exists between fiows within a 
fiux surface, which are first order in the pe/^i. 1 expansion, and radial transport fiuxes, 
which are second order in this expansion. This scale separation at least partially breaks 
down in the pedestal, and radial transport fiuxes compete with fiux surface fiows, even 
within a purely neoclassical framework. Our work includes this important effect, which 
strongly impacts the resulting flux surface flows. 

It is harder experimentally to measure the local bootstrap current density than to mea- 
sure the plasma flow. Since the current is just the difference in ion and electron flows, 
validation of neoclassical flow calculations would give confldence in bootstrap current pre- 
dictions. Impurity and main-ion flows have been measured and compared with neoclassical 
predictions in several experiments, with mixed results^"—. Neoclassical theory makes an 
absolute prediction for the poloidal flow but not the toroidal or parallel flows, since the 
latter are a function of d^o/dtp, and this radial electric fleld cannot be determined within 
the lowest-order axisymmetric theory^^. (Here, 27rilj is the poloidal flux.) The poloidal flows 
are largest in the steep-gradient transport barrier regions, yet these are precisely the regions 
in which the theory breaks down. For this reason as well, an improved nonlocal calculation 
of flows is sought to compare with measurements. 

Even in the limit of small collisionality, the form of the collision operator is crucial for 
determining the neoclassical flows, fluxes, and current. The collision operator rigorously 
derived from flrst principles is the Fokker-Planck-Landau operator—. In much analytic and 
numerical work, however, simpler "model" collision operators are used insteadi'^'^^"— . Model 
operators generally yield somewhat different results for all neoclassical quantities^ii^^, so in 
the following work the exact linearized Fokker-Planck-Landau operator is used. At the same 
time, it should be remembered that even the exact Fokker-Planck operator is only correct 
to 0(1/ In A) where In A is the Coulomb logarithm. 

Calculations of neoclassical quantities at realistic aspect ratio and with a realistic treat- 
ment of collisions require a numerical treatment. Local neoclassical computations with 
complete Fokker-Planck-Landau collisions were described by Sauter et al in Refs. 1231-125 



later extended to stellarator geometry in the code NEO described in Refs. 



26 



and 



and 



27. More 



recently, Fokker-Planck-Landau collisions have been implemented in other codes^^, including 



28. All of these 



a second code called NEO^ (unrelated to Refs. l26| andl27|), and in Ref. 
codes are radially local. 

In recent years, a number of numerical efforts have been undertaken to compute nonlocal 
neoclassical effects in transport barriers. Most of these efforts have used the particle-in- 
cell (PIC) approacb22i"— . PIC and continuum codes have differing treatments of collisions 
and boundary conditions, and face different numerical resolution challenges, so it is good 
practice to develop both approaches to verify they yield the same physical results. Some 
investigations of neoclassical effects have been begun in global continuum codes22i"— , but 
these codes use approximate collision models and are ultimately designed for turbulence 
studies, and very different algorithms have been used than the ones we use here. 

In this work, we present a new approach to computing global neoclassical effects. A 
continuum (Eulerian) framework is used, including the exact linearized Fokker-Planck- 
Landau collision operator. Our approach includes a general prescription for extending a 
local neoclassical code to incorporate nonlocal effects in a numerically efficient manner, by 
making such a local calculation the inner step of an iteration loop. Several terms related 
to the electric field must first be added to the local code, since in the pedestal these terms 
cannot be neglected. 

Our approach is not completely general, for while we allow the ion density scale length 
r„ to be ~ pe, we require the ion temperature scale length be > pe- The electron 
temperature scale length rx^ may be either ~ or > p^. The ordering r^. > r„ is satisfied 
in the pedestal on many present tokamaks, including DIII-D, JET, ASDEX-U, NSTX, and 
MAST—"—, when inter-ELM or non-ELMy (without RMP) profiles are carefully examined, 
though exceptions exist such as I-mode and EDA H-Mode in Alcator C-Mod^. Both entropy 
considerations^ and datai^i^ suggest rj; resists becoming as small as pe when collisionality 
is low, while r„ and r^^ are not similarly constrained. This suggests that future experiments 
with higher pedestal temperatures and lower coUisionalities may increasingly satisfy r^i > r^. 
The entropy argument is based on a more careful version of the asymptotic analysis above, 
showing that while r^. ~ pg would require j\ to depart strongly from a Maxwellian flux 
function, the same need not be true if r-r^ ~ pg or r„ ~ p^. Collisionless orbits radially 
average the ion temperature within a poloidal gyroradius, preventing strong ion temperature 
variation, while the density is not similarly averaged due to electrostatic confinement. The 



strong temperature gradient case for general collisionality requires use of the full nonlinear 
Fokker-Planck-Landau collision operator, giving rise to a kinetic equation that is nonlinear 
in /i. In the weak- 7]' case we consider, we will show it is still appropriate to expand about 
a Maxwellian flux function /mi and use the linearized collision operator. The E x B-driit 
nonlinearity associated with the poloidal electric field also becomes negligible, making the 
coUisionless part of the kinetic equation linear in 5/ = /i — fui- The full-/ strong- 7]' case 
will be considered in future work. Any more general full-/ nonlinear code must be able to 
accurately reproduce the weak-7]' limit, and since expansion about a Maxwellian is useful 
both numerically and analytically, it is worth understanding this limit in detail. 

Another simplification in this work is that we assume Bg ^ B, separating pe from the 
gyroradius p = Vi/Q scale. Without this approximation, the desired ordering r±_ ~ po would 
then imply the equilibrium varies on the p scale, so a drift-kinetic description would not be 
possible. This is not a serious limitation and is well-satisfied for the edge region, which is 
characterized by large safety factors. 

The primary finding in this work is that the ion flow is significantly altered in magnitude 
and direction relative to the prediction of local theory, and in particular, the flow's poloidal 
variation is qualitatively different. The poloidal variation of the flow is effectively determined 
by the requirement that the total flow be divergence- free. In conventional theory, the total 
flow is approximately given by a sum of parallel, diamagnetic, and leading-order E x B 
components, implying the poloidal flow must vary on a flux surface as Bq. However, in 
the pedestal, two other contributions to the flow divergence grow to become leading-order 
terms: the E x B flow of the poloidally- varying part of the density, and the radial variation 
in particle flux. As a result, the coefficients that multiply the ion temperature gradient in 
the parallel and poloidal flow are no longer equal, the poloidal flow no longer varies as Bq, 
and the flow may change magnitude and sign relative to the local prediction. These effects 
may be important to consider in any comparison between experimental flow measurements 
in the pedestal and neoclassical theoryii^i^. We present the details of one calculation of 
these effects at experimentally relevant aspect ratio and collisionality, considering a single 
ion species. 

In the following section we review the relevant aspects of local neoclassical theory. At 
the core of our global solver is a local solver, so in the next section we discuss in detail the 
local solver used. New comparisons to reduced analytic models are presented. Section [IV] 
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then discusses a Sf formulation for the global neoclassical problem in a transport barrier 
with a strong radial density gradient. Changes to the structure of the flow are discussed 
in section IVj and changes to the Pfirsch-Schliiter and bootstrap currents are calculated in 
section I VII Even in the 6f formulation, the kinetic equation is challenging to solve by direct 
numerical methods, so section IVIII introduces the operator-splitting initial-value-problem 
approach which reduces the dimension of the numerical problem to solve. Section IVIIII 
discusses the need for a sink term in the model and describes the sinks used. Results are 
presented in section HXf and we conclude in section |Xl 



II. DEFINITIONS AND LOCAL THEORY 

In the local case, the ion distribution function /; is approximately a Maxwellian with con- 
stant density ni{ilj) and temperature Ti(ip) on each flux surface: /mi = ni [m;/ (27rTi)]'^^^ exp (—miv'^ / [2Ti]) . 
To next order, fi = fui — Ze^ifui/Ti + /i where ^i{ip, 9) = ^ — $0) ^ is the electrostatic 
potential, $o(V^) is the flux surface average of $, and it can be shown |<I>i| ^ |$o|- The 
distribution /i is found by solving the following drift-kinetic equation: 

V||/i + K ■ Vij)dfMi/dij = (1) 

Here, ■ Vi/j = + v\/2) / {ytB'^)B x V-B ■ ViIj is the radial magnetic drift, and Ci is 
the linearized ion-ion Fokker-Planck-Landau collision operator. The derivatives in ([T]) are 
performed at fixed /i = miv\/ {2B) and total energy Wq = miv"^ /2 + Ze^o, so 

dfui 
dijj 

where x = v/vi. We ignore the 0{\/rnJmi) correction introduced by ion-electron collisions. 

It is sometimes convenient to apply the identity ■ Vip = {Iv\\/Q)V \\{v\\/ B) , where / 
equals the toroidal field i?^ times the major radius R, to rewrite ([1]) as 

vii^iig = Q{g + F} = Q{g} + Q {F} . (3) 

Here, F = —{Iv\\/il)dfMi/dijj, and 

9 = fi-F = f-fMi + ^e^i/Mi/Ti - F (4) 

Only a gradient can drive g, not gradients in or $o- This result follows from Ci{f ||/Mi} = 
0, so the dpi/dip and d^o/dip terms in ([2]) disappear entirely from Ci{F} and from ([3]). 
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Idpi ^ Ze d% ^ / 2 _ 5\ 1 dTi 
Pi dijj Ti difj \ 2 J Ti dijj 



/Mi (2) 



Once /i is found, the two moments of greatest interest are the radial heat flux 

■ V,) ^ (/ /, (!!|! _ . V,) ^ (5) 
and the parallel flow 

Here, fcg and k\\ are dimensionless coefficients defined by the right equalities in ([5])- ([6]), e is 
the inverse aspect ratio, and brackets denote a flux surface average: 

(A) = {V')-' / dOA/B-Ve (7) 
Jo 

for any quantity A where 

/27r I* 
de/B ■\/e = j die/ Be, (8) 

die is the poloidal length element, V = dV/dip, and 27rV{ip) is the volume enclosed by 
the flux surface. Also, = 4A/27rZ'^e^ni In A/ (^3^/rniT^^^^ = \f2v\ is the ion-ion collision 
frequency with vi the Braginskii ion collision frequency. The definition for in ([5]) turns 
out to be convenient as kq then has a finite limit as e — > and collisionality — 0. 

It can be shown using the following argument that the parallel flow must have the 
form ([6]) with k\\ constant on a flux surface. First, apply the operation J d^vi^ ■ ) = 
2TxBm^^^^o dv JJ"'^ ^^'^^^ d^{v/v\\){ ■ ) to ([3]). (Here, a = sgn(t;||).) This operation 
annihilates the linearized collision operator terms by particle conservation. Pulling V|| in 
front of the velocity integrals, the boundary term from the upper limit of the dfi integral 
vanishes in the a sum, leaving i?V||(/ d^vv\\g/ B) = 0, and so / d^vv\\g = n^XB where X is 
constant on a flux surface. Then applying n^^ f d^vv\\{ ■ ) to 1^, and noting the last term 
in ([2]) vanishes in the v integral, the flow must have the form 

cl f 1 dpi d^Q 



Vn = - + Ze — - + XB. (9) 

" ZeB \n;dip di) ) ^ ' 

Recalling g oc dTi/dip, and normalizing X by convenient constants, then the form results 
with k\\ constant on a flux surface. The form can also be understood from a fluid 
perspective, as follows. First, the leading-order perpendicular flow is V± = cB~^{d^o/dip + 
[Zeni\~^dpi/dip)B x Vip. Together with the mass continuity relation V ■ {uiV) = + 
0{niVipg/r'j_) and B x Vtjj = IB — R^B^VC for toroidal angle ( (which follows from B = 



VC X V'?/' + IV C), this implies (Q. The constancy of k\\ on a flux surface may be used as a 
test for any numerical scheme. 

The constant k\\ may also be understood as the magnitude of the poloidal flow Vq = 
V-eg = V^-ee+V\\B-ee/B = k\\cIBe{Ze {B^)y^ dT^/dtP where = (VCx V^/')/|VCx V^AI- 
This result arises because the Ex B and diamagnetic perpendicular flows cancel the dpi/difj 
and d^o/dip terms in ([6]) when the poloidal component is formed, leaving only dTjdip to 
drive poloidal flow. The coefficient k\\ arises again in the dTi/dip contribution to the parallel 
current, as will be shown in Section IVTl 



(j||fi> =(T„eo (^11 5) -Cjpe 



1 f dpe dpi\ C32dTc Cuk\\dTi 
^31 — I -r— + -1— I + 



Pc \d'ip dtp J Te dip ZTc dip 



(10) 



Here, cXneo, -^^si, -^^325 and £34 (defined in Section IVTl) are coefficients determined by the 
magnetic geometry and electron collisionality, affected by the ions only through the ion 
charge Z. 

The linearized Fokker-Planck-Landau operator for ion-ion collisions, needed to solve ([1]) 
or ([3]), may be written 



30F d 

dx 



xe 



■*(x^ 



_d g_ 

dx 



+ 3e- 
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H x^ d^G 

+ 



where z/d = (3v/vr/4) [erf(x) - ^(x)] /x^ ^ 
27r"^/2 e~^^dy is the error function. 



2nvf ' 2-711!^ dx"^ 
erf (x) — 27r~^/2xe~^^ /{2x^), erf(x) 



111 



d 



(12) 



is the Lorentz operator, and ^ = v\\/v. Also, H and G are the non-Maxwellian corrections 
to the Rosenbluth potentials, defined by V%H = —4:7rg and V^G = 2H, with the velocity- 
space Laplacian V|; = ir^ [{d/dx)x'^{d/dx) + 2L]. The last three terms in (ITT!) (those 
following 3e~^^) together form the "field part" of the operator. While historically this part 
of the operator is often replaced with ad-hoc models, here we retain the exact field terms. 



The concise form ( |TT|) of the field operator is derived in Eq. (7) of Ref. 
equivalent to the full linearized Fokker-Planck-Landau field operator. 



54j . and it is exactly 



III. LOCAL SOLVER 

The basic approach to solving the kinetic equation (jl]) with the full field operator is to 
treat H and G as unknown fields along with the distribution function g, and to solve a block 
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linear system for three simultaneous equations: the kinetic equation and the two Poisson 



equations that define the potentials. Figure [T] illustrates the structure o 
The approach is simi 



independently. Ref. 



ar to the innovative method described in Ref. 



28 



28 



this linear system, 
but was developed 



(a radially local code) is focused on the banana regime in which 
g = g{fi, v) is a function of two phase-space variables, whereas in the analysis here we wish to 
keep the coUisionality general, which means g depends on three or four phase-space variables 
(in the local and global cases respectively.) 

We may solve either ([1]) for /i or ([3]) for g. The operator and matrix are the same 
for the two approaches, but the right-hand side vector (the inhomogeneity) is different. 
The equivalence of the distribution functions obtained by the two approaches is another 
useful test of convergence. For the second approach, the inhomogeneous term in may be 
evaluated explicitly: 

- [^°^^"^' ^ '"^'^ ('^^ - ■ 

Deriving this result amounts to evaluating Ci{t'||f ^/mi}, which is done in e.g. Eqn. (C19) of 



Ref. 
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We discretize the Rosenbluth potentials by retaining a finite number of Legendre poly- 
nomial modes Pt{d- There are several motivations for this choice. First, the Legendre 
amplitudes of H and G fall off rapidly with £ since ~ Therefore only 2-4 modes are 
sufficient for convergence, although the code allows for the retention of an arbitrary number 
of modes. Secondly, the Legendre representation allows a convenient and efficient treatment 
of the boundary at large f , which can be understood as follows. The distribution function will 
be within machine precision of zero for v > 6vi, so it is wasteful to store g for this v region. 
However, H and G scale as powers of v rather than as e-(^/''-)', so they remain nonnegligible 
even for v > Qvi. (In fact, for general g, G increases with v.) However, with a Legendre 
representation H = Hi{v)Pi{^), we may exploit the fact that for v > fMax = 4 — 6 fi, 

the defining equation for H becomes VlH = 0, and so Hi = Aiv'^^^^^ + B^v^ . The physical 
solutions have = 0, and so the Robin boundary condition vdH^/dv + {i + l)i?£ = 
may be applied at ^Max to ensure oc v~^^^^\ In the case of V^G = 2if, there are 
four linearly independent solutions for G. Two are homogeneous solutions to V^G = as 
for H above, and two are particular solutions, which vary as v"^ times the homogeneous 
solutions. Thus G = J2T=oGi{v)Pi{0 where Ge = Gev-^^+^^ + D^v^ + Eev^-^ + Ft,v^+'^. 



Local DKE: 

M,V,i.-C{i;} 
= C{F} + S 

S\H + Ans =0 
V;G-2// =0 



J 



H 



C{F} 
+ S 



Vector of unknowns 



FIG. 1. (Color online) Block structure of the linear system for the local solver. A few rows are 
also reserved for boundary conditions. 

The physical solutions have Di = Fe = (see e.g. (45) of [14]), leaving one homogeneous 
and one particular solution. To accommodate both solutions requires a second order equa- 
tion as a boundary condition. Writing d'^Ge/dv'^ + avdG^/dv + bGe = 0, and inserting 
Gi oc v~^^~^^\ then Gi oc v^~^, yields two equations for (a, b), giving the boundary condition 
d^Ge/dv^ + {2i + l)v dGt/dv + {f - l)G^ = 0. 

The other boundary conditions applied are as follows: Hi = and = at f = for 
i > 0, dHi/dv = and dGn/dv = at = for £ = 0, = at = t'Max, dg/d^ = at 
V = 0, and dg/dv = at (v,^) = (0, 0). No boundary conditions are applied to g at C, = ±1 
(i.e. the kinetic equation is applied there with one-sided derivatives.) 

While it seems essential to represent the pitch-angle dependence of the potentials using 
Legendre polynomials, the distribution function itself need not be discretized in the same 
way, and there are many options available for the other coordinates, so a range of different 
discretization schemes were investigated. A choice of piecewise Chebyshev spectral coloca- 
tion and finite-difference methods of various orders were implemented for both the x and ^ 
grids. The spectral colocation approach is highly accurate for given grid resolution. How- 
ever, as the matrix is denser in the associated coordinate for these approaches, the solver 
slows more rapidly compared to finite differencing as the grid resolution increases. Thus, for 
satisfactory numerical convergence, high-order finite-difference methods are often preferable 
in practice. For discretization in 6, finite-difference methods of various orders and spectral 
colocation as well as a sine/cosine modal representation have been implemented. The modal 
approach is extremely efficient for the simple concentric circular flux surface model, in which 
case the matrix is sparse in 6. However, for shaped geometry, the matrix becomes dense in 9 
for the modal approach, so the colocation approach is both more convenient and similarly ac- 
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curate. In shaped geometry, despite the accuracy of the spectral approaches, finite-difference 
differentiation again typically gives satisfactory convergence in less time due to the sparsity 
of the matrix. 

The linear system may be solved using a sparse direct algorithm; a Krylov-space iterative 
solver may be much faster, but convergence of the algorithm then requires an effective 
preconditioner. One successful preconditioner is obtained by eliminating the off-diagonal 
blocks in Fig. [1] as well as the off-diagonal-in-x terms in the energy scattering operator and 
boundary conditions. If high-order finite difference derivatives are used in x, convergence 
typically also requires that a constant ~ ua be added to the diagonal of the kinetic equation. 
We find the generalized minimum residual method (GMRES) does not converge consistently, 
while the stabilized biconjugate gradient and transpose-free quasi-minimal residual methods 
are more reliable. 

Several issues regarding null solutions and symmetry properties of the distribution func- 
tion are discussed in Appendix Rl 

Figures [2] and [3] show typical results of the local code, plotting the fiow and thermal 
conductivity coefficients k\\ and kg as functions of aspect ratio and collisionality. Although 
the code can use general shaped geometry, for all plots in this paper we use the standard 
concentric circular fiux surface model B oc 1/(1 + e cos 6) and b ■ V6 = constant, to facilitate 
comparison with previous literature on neoclassical theory. We may then take the definition 
of the ion collisionality to be z/^, = i'ii/{e^^'^v{V\\6). 

Several approximate analytic formulae are also plotted. The Chang-Hinton formula for 
the heat fiux— 



kn 



0.66 + 1.886^/2 _ 2 

^72 



5^ 



1 + 1.03Z/*' +0.31Z/* 
and the formula of Sauter et al^ 

1 ^ -1.17/c 



1 

52 



0.58z/,e 



0.74z/,e3/2 



1 

52 



1 



(14) 



k\\ 



1 + O.Sy^ VI - 0.22/t - 0.19/t^ 



+ 0.25(1 - /^)v/];: ) +0.315z/,Vf 



1 



0.15z/2jt' 
(15) 



apply to arbitrary aspect ratio, plasma shaping, and collisionality. (Note a in Ref. 
—k\\ in our notation.) Taguchi's formula for the heat fiux— 



25 



fc 



/c + 0.462/t. 



equals 



(16) 
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FIG. 2. (Color online) The flow coefficient k^^ defined in ([6]) for concentric circular flux surface 
geometry, as computed by the local code (solid) and the Sauter et al formula psp (dashed). Squares 
show the Helander-Sigmar formula (|17p . an approximate analytical treatment of the z^* — t- limit, 
for the same four values of e. The arrow indicates the known analytic z^* — >■ 0, e — )• limit 1.17. 



and a formula for the flow coefficient derived on p.216 of Ref 



fc|| = 1.17/c/(/c + 0.462/t: 



■.a 



(17) 



are applicable at arbitrary aspect ratio and shaping in the limit of small i^*. An expression 
equivalent to the latter result was also given previously in Eq. (28) of Ref. Q. (There, 
a2 = 0.6562A/e and ai = 1.173. Taking ft = lA6^/e for circular flux surfaces yields the 
same result as ( 1T7|) .) Here, /t = 1 — /c and 

"1/ -Bma 



fc 



(18) 



^0 VI - A5 

As Figure [2] shows, ( fT7|) does a reasonable job of predicting the low-collisionality limit of k\\. 
The analytic result k\\ = 1.17 obtained using a momentum-conserving pitch-angle scattering 
model collision operator is indeed the limiting value for z/* — and e — t- as expected, 
but e must be < 0.01 for this value to be a good approximation. Figure [3] shows Taguchi's 
formula (fT6l) is extremely accurate. The Chang-Hinton formula (fT4l) is less accurate but it 
correctly captures the trends with e and u^. 

The local code was also compared with published results from the Fokker-Planck code of 



Ref. 



221 : the remarkable agreement of the codes is shown in Figure |H 
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FIG. 3. (Color online) The thermal conductivity coefficient kq defined in ([5]) for concentric circular 
flux surface geometry. Solid curves are computed by the local code. Dashed lines indicate the 
Chang-Hinton formula (|14|) for the same four values of e. Squares show the Taguchi formula (|16p . 
an approximate analytical treatment of the z/* — )■ limit, again for the same four e. The arrow 
indicates the known analytic z^^, — )• 0, e — )• limit 0.66. 




FIG. 4. (Color online) The local code agrees to high precision with the Fokker-Planck code of Ref. 



22l for both the flow (a) and heat flux (b) coefficients. Calculations are shown for e = 0.1. Data 



reproduced with permission. 
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FIG. 5. (Color online) a) Particle orbits, i.e contours of magnetic moment /x, for e = 0.3. b) The 
"collisional response" distribution g at v = Vi for i'^, = 0.01, showing ~ for the trapped region 
as predicted by banana-regime analytic theory, c) A slice of g at 9 = shows the boundary layer, 
which is narrower at larger v due to the lower effective collisionality. 

Another set of transport coefficients arise in the analysis of the electrons. The radial 
particle and electron heat diffusivity are ~ ^Jm^jm^ smaller than the ion heat transport 
and are always dominated by turbulent transport in practice, so we will not discuss them 
further. Of greater interest are the electrical conductivity and bootstrap current; these 
quantities are discussed in Section IVll 

The distribution function obtained by the local solver has several noteworthy features. 
In the z/* <^ 1 limit, analysis shows the g piece of the distribution function should vanish in 
the trapped region of phase space^i^, and a boundary layer exists between the trapped and 
passing regions^. These properties are reproduced in the code, as illustrated in Figure O 
The thickness of this boundary layer increases with collisionality. Although the collisionality 
z/* is typically defined in terms of the thermal speed fi, each value of v in phase space 
effectively has its own collisionality given by vj^j {ve'l''^h ■ V6'), with faster particles being 
less collisional. As shown in the figure, the boundary layer indeed grows narrower with v. 
In Figure Oc, g at each v is scaled to go to 1 at ^ = 1 for clarity. Also notice in Figure |5la-b 
that g is nearly constant along particle trajectories, as it should be. 

IV. GLOBAL KINETIC EQUATION 

Under what circumstances does the ion distribution remain close to a Maxwellian flux 
function, i.e., is the approximation /i ^ /Mi(^/') of section HTl valid? The magnitude of the 
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correction to the flux-function Maxwellian may be estimated from the local theory roughly as 
f^r^Fr^ {pe / ^ 1.) fui where r±_ is the scale length of variation in density and/or temperature. 
In a pedestal, since pe/fi, ~ 1, then /i ~ /mi and the neoclassical expansion breaks down. 

However, a more careful estimate reveals there is a regime in which the near-Maxwellian 
assumption is still valid^. To deflne this regime, flrst write fyn in an equivalent form: 

3/2 

/Mi = 



-P I ) (19) 



where 

r7(V') = n,{i,) exp (Ze$o(^)/Ti(^)) (20) 

and again Wq = miv"^ /2 + Ze$o is the leading-order total energy. Then the derivative 
dfui/dip (at flxed Wq) that determines the magnitude of /i is 

dfui _\ldr] f 3\ 1 rfTil 

U# + r°"2j7]#J^^^ ^^^^ 

(equivalent to ([2]).) In this form, it is apparent that the magnitude of dfui/dtp is determined 
by tt and r^, the scale-lengths of and rj, but not directly by r„, the scale-length of density. 
Therefore fi/fui may be small compared to unity even when r„ ~ pg as long as tt and 
are > p^- 

This "weak-TJ' pedestal" regime is the ordering^"^ we shall consider for the rest of the 
analysis: 5^1 where 6 = pe/rx, Pe/^rj ~ S, and pg/rn ~ 1. This regime is useful for two 
reasons: the collision operator may be linearized, and, as we will show, the poloidal elec- 
tric fleld may be decoupled from the kinetic equation, eliminating the E x B nonlinearity. 
Therefore the kinetic equation is linear in /i. For tt ~ pe and/or pg, a full-/ nonlinear 
kinetic equation must be solved, retaining both the coUisional and E x B nonlinearities. 
Notice <^ 1 implies the ions are electrostatically conflned {d^o/dip ~ —[Zeni\~^dpi/di/:) 
with {Ze/Ti)d^o/dip ~ 1/ (RBgpg), so Ze^o/Ti ~ 1. Due to this ordering for the electric 
fleld, the Ve ■ V/ term in the kinetic equation, neglected in conventional neoclassical calcu- 
lations, becomes comparable to the t>||V||/i streaming term. Therefore, although /i <C fui 
in the weak-7]' pedestal, conventional neoclassical results must still be modifled. 

Now, consider the full-/ drift-kinetic equation^ 

t'||V||/ + Vd- V/ = C„i{/} + 5 (22) 

where 5" represents any sources/sinks and C^i is the nonlinear Fokker- Planck- Landau opera- 
tor. As pointed out by Hazeltine^, fl22|) may be derived recursively, and so its validity does 
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not require |u||V||/i| ^ \vd ■ V/i|. Since fi ~ fui to leading order, Cni niay be approximated 
with Ci, the operator hnearized about fui- For the drift velocity it will be convenient to 
use = (f||/fi)V X (v^^b) (discussed in appendix [B|) where the gradient acts at fixed /i and 
W. We make the ansatz $i ~ 5 $o and d^i/dip ~ 6 d^o/dtp, and we will show shortly that 
these assumptions are self-consistent. We define /i by /i = fui — (^e$i/Ti)/Mi + /i, and 
change the independent variable from W to Wq- Neglecting several $i terms that are small 
in S, 

{mb + v^) ■ V/i - = -^^d ■ + S. (23) 

The contribution from $i to ■ V6 is 0{6) smaller than the $o contribution, and {ve ■ 
Vip)/{vm ■ V?/') ~ Ze$i/Ti ~ 6 where = — ve is the magnetic drift, so $1 may be 
entirely neglected in and (!23|) . We therefore approximate in (!23|) with Djo = i'm + 'i'_Eo 
where Veo = cB~^B x V$o- To evaluate $1 we may use the adiabatic electron density 
response ng + (e$i/Tc)ne, where nc{ip) = Zni{ip) is the leading-order electron density, with 
quasineutrality to obtain 

e<l>i/Ti = (Ti/Te + Zy' J f,. (24) 

Hence, as /i ~ Sfui, the ordering ansatz for $1 above is self-consistent. As both the 
coUisional and E x B nonlinearities are thus formally negligible, ( 123|) is completely linear. 

Just as in the local case, it is convenient to define the collisional response part of the 
distribution function g using (j4]). Eliminating /i in ( l23l) in favor of g, a pair of terms 
cancels. We also drop the resulting i>do' V [{Iv\\/Q)dfMi/dip~\ term because -Udo' V(/f||/-B) oc 
V X [{vii/B)B] ■ V{Iv\\/ B) = exactly and Vdo • "^{dfui/dip) is small in 6. Thus, we obtain 

{viib + -.^do) ■ - Q{g} = Q{F} + S. (25) 

The advantage of this second form is that it makes clear that gradients in ni, $0? and/or 77 
cannot affect the g part of the distribution function - only a T[ gradient can drive g. The 
logic is the same as in the local case: Ci{f||/Mi} = 0, so the only gradient surviving in the 
inhomogeneous drive term (fT3l) is dTi/dip. This property is obscured in the form (l23l) . While 
the independence of g from dni/dip and d^o/dip was well known previously for the local case, 
it is noteworthy that this property persists in the weak-T/ pedestal case considered here^. 
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V. CHANGES TO FLOW STRUCTURE 



Two noteworthy differences between the local and global analyses are that the parallel 
ffow coefficient k\\, as defined in ([H]), no longer needs to be constant on a flux surface, and it 
no longer also describes the poloidal flow. To see the first of these points, we may apply the 
operation J d^v to the kinetic equation (!25l) . as detailed in appendix [Bl The resulting mass 
conservation equation (ignoring sources) is 

V • (^J dh{vilb + VEO + Vm)^^ = 0. (26) 

Recall from section [III that in the local case, the v\\ term in fl2Bl) dominates the others, 
implying J d^vv\\g oc B. This result was crucial for proving the constancy of k\\ on a fiux 
surface, for the dT^/dip term in the parallel fiow is precisely j d^vv\\g. However, in the global 
case, fl26l) indicates that j d^vv\\g need not vary on a fiux surface in proportion to B, so the 
proof for the constancy of k\\ breaks down. 

These same results can be derived from a fiuid perspective, making no reference to the 
drift-kinetic equation, starting instead from the fiuid mass fiow 

T = j d\> {v\\h + VE + v„,) fi-V xM (27) 

where M = b j d^vfiv\/Q. Equivalent ly, 

r = r||b + r^ + rdia, (28) 

where Te = J d^v fiVE, Tdia = c{ZeB'^)^^B x V is the diamagnetic fiow, 

p±{ I — bb) + p\\bb, p± = rriij d^v fivj_/2, and p\\ = m\ j d^v fiv'^y The equivalence of fl271l 

and (128|) follows from 

j d'v fiV^ - V X M = Tdia, (29) 

which may be derived^ using the more accurate drift = v'^^Q^'^b x k + vj_(2QB)^^b x 
Vi? + v'j_{2QB)^^bb ■ V X b. (This drift is identical to our earlier expression to lead- 
ing order in (3 <^ 1.) Notice b- (l25l) with (jl]) and Vji = F ■ b/rii gives with = 
Ze {B"^) [clriiB dTi/dip)~^ j d^vv\\v as before. 

We now impose mass conservation V-F = 0, substituting (jlj) into ( 1271) . applying BxVip = 
IB - R'^B'^VC, and noting / d^v v^fui = c{ZeB'^Y\dpJ di))B x + V x [cp-.b/ {ZeB)]. 
Cancellations occur to leave Ti + T2 + T3 + T4 = where Ti = V ■ / d^v{v\\b + veo + v^)g, 
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T2 = -V ■ / dMi^Eo + v^)Ze^jMi/Tu T3 = V ■ / cPwEifmi = V ■ {cmB-^B x V$i), 
and T4 = V ■ / d^v VEi{-ZeTr^ fui^^i + F + g) where vei = cB'^B x V$i. As Vri; = 
— (Zeni/Ti)V$o + 0(5) in our ordering, T2 and T3 cancel to leading order in 6. It can be 
verified that the terms in T4 are 0{6) smaller than the terms in Ti, so to leading order, Ti = 0, 
which is precisely f l2B]) . but re-derived from a fluid rather than drift-kinetic perspective. The 
fluid analysis thereby confirms k\\ is no longer constant on each flux surface. Compared 
to the fluid analysis in the conventional ordering, reviewed following (|9]), it can be seen 
that two new contributions to mass conservation become important: E x B convection 
of the poloidally varying density, and radial variation of the particle flux or (equivalently) 
diamagnetic flow. Even though \v eo\ Vi, E x B convection of the density carried by g 
matters for mass conservation fl2^ because the parallel flow only enters multiplied by the 
small factor Bg/B. And although ^ Pi I , the radial derivative in V ■ T means the 
next-order correction to the diamagnetic flow (or equivalently the radial neoclassical 

flux) must be retained to accurately compute V ■ T. 

The poloidal fluid flow Ve is found by computing Ve = T ■ eg/rii, using (1271) or (l28ll . 
Plugging (jlj) into efl-( l27l) . several cancelations occur, leaving 

Ve = 4^ [d^vv,g+'4^'^ [d^vg-^.Vx [ d\,$-gb + ^ ■ [ d\, v^g (30) 



Bui J B'^rii di/j J rii J n\ 



^ V- 



V2 V-i 



Ze$i cIBe 



d^o Ti drii 
dip Zerii dip 



SxV*,.e,-(-— + 



So far no terms have been dropped. We now order the terms using the orderings developed 
in Section HYI Using {Ze/T{)d^Q/dip ~ l/{RBepg) it can be verified that Vi ~ V2. It can 
also be verified that each term following V3 is 0{5) smaller than Vi ~ V2 using J d^v g ~ 6ni, 
Ze^i/Ti ~ 5, V<I>i ■ B X Be ^ -IBed^i/dip ~ 5IBed^o/dip, and noting the quantities in 
square brackets cancel to leading order. 

It remains to evaluate V3. The leading order contribution comes from the radial gradient 
of the integral of g, since only this derivative has the short scale length pe. Thus, we obtain 

riiB 

In the local case, only the v\\ term arose in the analogous integral for Ve- In the pedestal we 



i^ibf) ' ^ m r^i' ■ '''' 



may define a normalized poloidal flow 

ke = VeZe (B^) /{cIBe dT-Jdip) (32) 
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so ke k\\ in the local limit. The property Vq oc dTi/dip from conventional theory persists 
in the pedestal, due to fl3Tl) and g oc dTi/dip. 



VI. ELECTRON KINETICS AND PARALLEL CURRENT 



The orbit width for electrons is ~ ^J^m^ml thinner than that of the ions, so direct finite- 
orbit-width effects for electrons may be neglected. However, the electrons are affected by 
modifications to the main ion flow. To demonstrate this point, and to show applications of 
the local Fokker-Planck code to electron quantities, we now analyze the electron kinetics. 
Since the particle and electron heat transport are essentially always dominated by turbulent 
transport, we focus here instead on the neoclassical conductivity and bootstrap current. 
Though the analysis below uses the pedestal ordering, conventional results for the parallel 
current are exactly recovered in the appropriate limit of the expressions derived here. 

Using the gauge of Appendix O) the electron kinetic equation may be written 

{v\\h + + ve) ■ (V/e)^ - ev\\ {E\\B) B {B^y' dfjdw = a (33) 

where v^^ is the electron magnetic drift, w = — e$ is an independent variable, and 

Ce is the total electron collision operator. We assume /e ~ /mc where 



/Me = ne{lp) 



3/2 / 2 



Then Ce = Cee + Cci where C^e is equivalent to ( ITTl) but with ion quantities replaced by 
electron quantities, Cci{/ci} ~ UciL{fci} + fMe^eimcV\\Vi\\/Te, = 3^/7r /{ATdX^), x = v/v^, 
Ve = \f2Tjm'^, and Tci = 3y'm^Te^'^^/(4A/27rne2'e^ In A). We write /e = /mb exp(e$i/Tc) + 
mcV\\y\\]yiclTc + h and solve for h. We also make a change of independent variables in the 
kinetic equation to Wq = — e$o- Using ([6]), the leading terms in 6 and \fmyjm\ are 

W||V||/io + ^^me ■ V/Me + /mc (^^^ + ^"^^ '^'H^H ^^^^ " C'eo{/io} + Z^ei^{/io} (35) 

where ilg = —eB/{m(.c), ho is the first term in a series h = + hi + . . ., and the 
inductive term has been taken as higher order. The solution to ( 135|) may be written 
ho = cle~^{hpdp/dip + hT^UedTe/dip) where p = Pc + Pu and hp and hr^ are the solutions to 

D,hp = -fuen^x^l + e)B-^V\\B, (36) 
Dehn = -fMen^x^x' - 5/2)(l + e)B~^V\\B, (37) 
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with Dc = v^Vw — Ccc — J^ciL. Recalling e$i/Ti ~ S, the 0{S) terms in the kinetic equation 



are 



Dh - -^^^^11 V|| f ^^f^ ]+Vei- V/mc + e<^iV^, ■ V{fMe/%) 



where v^i = cB~'^B x V$i, and we have assumed <^ 6. The solution may 

be written hi = —hEe~^ (^E\\B^ — cIe~^nihT-^ dTi/ dip — pQcP{dnc/dil)){dTi/dil))h^/e where 
Po = ViTriic/ {ZeBs.^), and B^^ = (-B^). Here, He, fiT^, and are the solutions of 

DehE = fueeX"' {By' Bv\\, (39) 
= fueu^'x^ {BY' [^||(3e' - l)V||fi - 2eBV\\k\^ , (40) 

Note in the local case, V||/c|| = so the last term in PUI) vanishes and /it; oc k\\. The Z?e 
operator, which is radially local in that ip is merely a parameter, may be inverted numerically 
for the right-hand sides f l36|) - fl371) and fl39|) - PT|) just as described in Section |TTl] for the similar 
ion operator f || V|| — Ca. Then the parallel current is 

Jll = ZeuiViii -e j d\ fe = -e j d\^ v\\h. (42) 

Now consider the result of applying J d^v = a dv ^^'^^^ dfi Bv /v\\ to (!39|) . 

This operation annihilates both the right-hand side and the collision operators in De, leaving 
[d/dO) J d^vv\\hE/B = 0. Therefore the flow carried by /i^; is J d^vv\\hE = ueB for some 
flux function a^. The same logic applies to f l37|) . so j d^v v^hx^ = ax^B for some flux function 
ax^. Applying j d^v to ( 136|) . the right-hand side is not annihilated this time, and we instead 
find J d^vv\\hp = B^^ + apB for some flux function Op. Lastly applying j d^v to the first 
equation in ( HOl) and to ( HTl) . we obtain / d^vvwhr^ = ax^B + k\\B / {B"^) and / d^vv\\hii> = 
acs>B — ng/{ZB) where and are flux functions, Ug = Ti^pQlnidTj dip)~^ f d^v g is the 
0(1) normalized density carried by g, and we have invoked (Elj). Thus, the dT^/dip term in 
the parallel current varies poloidally oc 5 in the local case where k\\ is constant, but not in 
the global case where k\\ varies. 
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Putting the pieces together, the total parallel current is 

cl dp cIncBk\\dTi pQcPugdricdTi 

where a is another flux function. Multiplying this equation by B, fiux-surface averaging, 
and substituting the result back into fj43|) . we obtain 



c^dp_f^^ \ cln^B dTj f (B^hY 
'Bd^[jB^~ J^YlB^di^y "~ (52) 

Poc/^ due dTi / (ug) B 




Z dip dllj \ B J ' • (^^) 

The dp /dip term is the standard Pfirsch-Schliiter current, and the {i\\B^ term is the Ohmic 
and bootstrap contribution. However, the k\\ and Ug terms are new in the global case, 
vanishing in the local case where k\\ is constant and \Wnc\pe ^ 1- We may write the Ohmic 
and bootstrap contribution as 

/ . d\ d\ t f ^-ii dp Cs2dT^ Ct, dTi CnTpal dn^dTi\ 

{j,B) = .„eo {E,B) - dp. + — — - - ^^^^ ) (45) 

where (Xneo = {Bjd^vv\\hE), £31 = {B J d^v v\\hp) , C32 = {B J d^vv\\hT,), Lt, = 
i^B J d^vv\\hT-), and £„t = (^B J d^v v\\h^). The CnT term in ( 145|) is new, becoming 
negligible in the conventional case. For the local case of constant k\\, where Ct^ oc k\\, it is 
useful to define £34 = C-Tjk\\ so £34 is completely independent of all ioii quantities except Z. 

25l . Interestingly, the 



The definitions of (Tnoo, £31, £32, and £34 here are consistent with Ref. 
new Ug terms in ( 144|) and (143|) and the new CnT term in ( H5|) are quadratic in the gradients. 

Figure |6] shows these coefficients of the bootstrap current and the conductivity as calcu- 
lated by our code for the local limit A;||=constant, using the circular flux surface model and 
Z = 1. The conductivity has been normalized by the parallel Spitzer value. The analytic fits 
to numerical calculations of the coefficients by Sauter et al are plotted for comparisons^. The 
horizontal coordinate in these plots is i^^e = i^cd (e'^^^ \/ '^Tc/rajj ■ V9) , which is 1/a/2 smaller 



than the z/*(, defined in Ref. 
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We find the Sauter expressions give an excellent fit to the 
coefficients in the banana regime, though there is some discrepancy at higher coUisionality 
when e > 0.1, the same pattern observed in Figured The reason for the discrepancy is 
unclear, since the fundamental kinetic equations and collision operators we use to generate 
figures [2] and [6] are identical to those solved by CQLP, the code to which the Sauter ex- 
pressions are fit. (CQLP uses an adjoint method whereas our results do not, though this 
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difference should not affect the physical results.) We have verified the difference persists 



when D-shaped Miller equilibrium is used, and the code of Ref. 
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produces identical coef- 



ficients to ours. As shown in figure [71 the difference between our coefficients and those of 



Ref. 
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can lead us to predict a reduced total bootstrap current density in the pedestal for 
experimentally relevant plasma parameters when > 1. This difference is primarily due 
to our lower £31, which multiplies the large dp/ dip term. When z/*c < 1, our prediction for 
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the total bootstrap current density becomes indistinguishable from that of Ref. 

As with the flow, the total current vector j remains divergence-free in the pedestal or- 
dering: 

= V-j = V-(j||b + c5-2Sx V-tfs) (46) 

where tf . = mjd^vvv{U + /,) is the total diagonal anisotropic stress, including 0(p) 
and 0{5p) terms. Equation fjlB]) can be proved from (jUj) and fl2Bl) using (1211) . /e ~ /mb + 
^^ifuc/Tc, and (l29l) .) Equation (H6l) indicates the new k\\ and Ug terms in (jSj) arise for the 
same fundamental reason as the conventional Pfirsch-Schliiter current: a parallel current 
must flow to maintain V ■ j = given the perpendicular diamagnetic current. In the 
pedestal, the diamagnetic current associated with the poloidally varying pressure becomes 
large enough to modify the parallel current on the level of the dT^/dip terms. 



VII. GLOBAL NUMERICAL SCHEME 

It is equally valid and equally numerically challenging to solve either (l23ll or (l25l) . For 
the rest of the analysis here we discuss the case of ( !25|) for definiteness. 

As we are interested in a narrow radial domain around the pedestal, we assume J, -B, and 
V||6' are independent of ip for simplicity. These approximations are also convenient as they 
make I'm ■ = exactly for our form of the drifts. For simplicity, we also take rj and Tj to 
be constant over the simulation domain. The one place where dT^/dip must be retained is 
in the inhomogeneous term, since the drive is oc dTi/dip. As the kinetic equation is linear, g 
may be normalized by dT-Jdip, while every other appearance of Ti is treated as a constant. 

Both versions of the global kinetic equation fl25]) and fl2^ resemble their local counterparts 
dl]) and ([3]), but with the additional I'do' V term in the unknown. Due to the radial derivative 
in this term, the radial coordinate no longer enters the kinetic equation as a mere parameter, 
meaning the problem is now four-dimensional: g = g{ip, 6, n, Wq). In these original variables, 
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FIG. 6. (Color online) Parallel current coefficients defined in ([45]) computed by our local code 
(dots, connected by solid curves). Dashed curves show the semi-analytic formulae of Sauter et al^. 

23 



a) 1 
0.5 




b) 



10 

10° 




c) 

g 5 



d) 



r 




1/5 original colllisionality 



Sauter 
This work 




original colllisionality 



Sauter 
This worl< 



5x original colllisionality 



Sauter 
This worl< 



0.85 



0.9 



0.95 



FIG. 7. (Color online) a) Model profiles resembling the DIII-D measurements in figure (8) of Ref 
44 and a plausible q profile yield the collisionality profiles in b). c) The local code described here 



25l for lower collisionality. 



predicts identical total bootstrap current density to the formulae of Ref. 
d) We predict somewhat lower current density at the real collisionality, due primarily to the dis- 
crepancy in £31 shown in figure [Hb. e) At higher collisionalities, the differences become significant. 



the allowed range of each coordinate depends on the other coordinates in a complicated 
manner. For numerical work it is therefore convenient to change the independent variables 
from (/i, lio) to (^^,0 so the coordinate ranges become coordinate-independent. In these 
variables, the kinetic terms in ( 123|) and (!25|l become 

da 

{v\\b + t;do) ■ {Vg)^,w, = Ko{g} + Keig} + ■ V^-^ (47) 



where 



^o = ^^||(V||^)|-r;^i^(V||i?)| (48) 
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is the drift-kinetic operator implemented in conventional neoclassical codes^i^, and 

. . w| . «c.^ <yS(V„.)| - ... . ,4.) 

consists of new terms proportional to the radial electric field. In a pedestal, not only is the 
dg/dtp term in f H7|) important, but the terms in Ke also become equally important. The 
aforementioned ordering for d^o/dip implies each term in Ke comparable in magnitude to 
Kq. Physically, the latter two terms in (H9|) are essential for maintaining conservation of n 
and total energy as a particle's kinetic energy changes during an orbit. This kinetic energy 
changes because the electrostatic potential seen by the particle varies over an orbit width. 

We choose ni{ip), which determines $0 = {Ze)~'^T[ln{ri/ni). On either end of the radial 
domain, we take ni{ip) and ^o{ip) to be uniform for a distance of several pe, as illustrated 
in figure Oa-b. In this way, the distribution function will approximate the local neoclassical 
solution at the radial boundaries, so local solutions can be used there as inhomogeneous 
Dirichlet boundary conditions. 

As the kinetic equation ( l25l) is linear, it may in principle be solved numerically using a 
single matrix inversion. Indeed, this is the approach traditionally adopted by local neoclas- 
sical codes^"— including the one described in Section lllli However, this approach is al- 
ready somewhat numerically challenging for the local problem due to the three-dimensional 
phase space, as the matrix has dimension {N6 Nv NC,) x [N6 Nv N^), where NO, Nv, 
and A^^ are the number of modes or grid points in the respective coordinates. In the 
nonlocal case, the additional spatial dimension means the matrix size must increase to 
{Nip NO Nv N$^) X [Nip NO Nv NC,) for Nip radial grid points, making such an approach 
much more time- and memory-intensive. Therefore we seek an alternative method. 

In the new approach proposed here, a derivative with respect to a fictitious time dg/dt 
is first added to the left-hand side of (l25l) . For reasonable initial conditions and boundary 
conditions, g should evolve towards an equilibrium since the equation (125|) is dissipative. 
However, an explicit time-advance requires very small time steps for stability due to the 
many derivatives in the kinetic equation, and an implicit time-advance would require the 
inversion of a matrix just as large as for a direct solution of the original time-independent 
equation. 

An effective solution is to employ the following operator-splitting technique. Consider 
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the following series of two backwards-Euler time steps: 

'""7,~" +/fNL{ft.a/2,} = o. (50) 

7;""^" + = QIF] + S. (51) 

where Ki^ = Kq + Ke — Ci is the "local operator" and /Tnl = iv^-'Vip)d/dip is the "nonlocal 
operator." When summed together, gt+(i/2) cancels, leaving an equation that is equivalent 
to first order in At to a backwards-Euler time step with the complete operator /Tnl + -^l- 
However, each of the steps fl50|) - fl5T]) are much easier than a step with the total operator 
because the dimensionality is reduced: e.g ip is only a parameter in (150|) . so this step requires 
the inversion of Ntp matrices, each of size {N6 Nv N^) x (A^^ Nv N$,). Also notice that the 
local operator at each radial grid point need only be Lt/-factorized once, with the L and U 
factors reused at each time step for rapid implicit solves. The same is true of the nonlocal 
operator at each v and ^. 

Several higher-order operator splitting schemes were explored, but none were found to be 
stable for the equation here. 

The procedure outlined here provides a general recipe for extending a conventional neo- 
classical code into a pedestal code. A conventional neoclassical code inverts an operator 
Kq — Ci, i.e. many of the terms in Ki^, so minor modifications would allow such a code to 
carry out the local part of the time advance. The modifications necessary are adding the 
electric field terms Ke and adding the diagonal associated with the time derivative. The 
resulting operator is then iterated with the nonlocal operator. 

For the results shown here we employ a piecewise-Chebyshev grid in ip with spectral colo- 
cation differentiation. A tiny artificial viscosity is required at the endpoints for numerical 
stability; the magnitude of this viscosity may be varied by many orders of magnitude with 
no perceptible change to the results. Inhomogeneous Dirichlet radial boundary conditions 
are imposed, with the distribution function at these points taken from the local code. For 
completeness, we have also tried upwinded high-order finite differences for radial differenti- 
ation, with the upwinding direction opposite above and below the midplane, corresponding 
to whether drift trajectories in the region move towards increasing or decreasing ip. For our 
sign convention, the magnetic drifts are downward, so the inhomogeneous Dirichlet radial 
boundary condition must be specified above the midplane at large minor radius and be- 
low the midplane at small minor radius. This radial discretization scheme gives equivalent 
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results to the Chebyshev method, but it requires more grid points for convergence, and a 
numerical instability tends to arise at large times. 

VIII. NEED FOR A SINK 

In order to reach equilibrium, it is essential to include a heat sink. This requirement may 
be understood physically as follows. As we take the scale-lengths at each radial boundary 
to be large compared to pe, the heat flux into the volume at small minor radius and the heat 
flux out of the volume at large minor radius are determined by the local neoclassical result 
These fluxes are different due to the different densities at the two boundaries, and so 
net heat will constantly leave (or enter) the simulation domain. More rigorously, as shown 
in appendix |B| the (J d^v{ ■ )) and Jjf™'"' c??/' V ( J d^v{miV^ /2){ ■ )) moments of the kinetic 
equation in steady state give 

(52) 

d^v gv^ ■ Vijj^ 
fi'^'^s)^ (53) 

The first equation represents local mass conservation, and the quantity following V' is the 
particle flux. The particle flux is exactly zero in the local limit, so it vanishes at the 
radial boundaries, and so in the absence of a source/sink, it must vanish everywhere in the 
domain. In the second equation, representing global energy conservation, the first term is 
the difference between the heat into and out of the domain, and the second term represents 
change in electrostatic energy associated with particle flux. If 5 = 0, then the latter two of 
the three terms in (153|) vanish, but the first term is nonzero because the heat fluxes at the 
two radial boundaries are generally unequal. This contradiction proves the kinetic equation 
has no steady-state solution without a sink S. 

In a real pedestal, there will be a divergence of the turbulent fluxes, which would act as 
a sink term in the long-wavelength (drift-kinetic) equation we simulate here. Determining 
the phase-space structure of this turbulent sink term from flrst principles is an extremely 
challenging task, beyond the scope of this work. We therefore use a variety of ad-hoc sink 

27 



1 d 
V 



dhiS 



V'{ I dhg^v^-ViP 



^min 



''/'max 

+ Ze I dipV'^ 



dtp V 



terms, and we find the simulation results are only mildly sensitive to the particular choice 
of sink. 

The standard sink we use is 

S = -1 {giO + 9{-0) (54) 

where 7 is a constant. The sum over signs of ensures that S vanishes exactly for an 
up-down symmetric magnetic field in the local limit, due to the parity of the local solution 
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for global 



discussed in appendix \M This sink is quite similar to the one described in Ref . 
6f gyrokinetic codes. The constant 7 may be varied by several orders of magnitude without 
major qualitative changes to the results. 
Another option we consider for the sink is 

S = -7m (^1) /Mi - 7p (Pi) - /Mi (55) 

where and 7p are constants, rii = J d^vg, and Pi = J d^v {miv'^ /3)g. The first term in 
( l55l) dissipates any mass in g, while the second term dissipates any energy in g. 



IX. RESULTS 

Figures IHIITO] show results of the global calculation for a pedestal with e = 0.3. The 
simulation domain consists of an annular region in space, i.e. an interval in ip. The density 
varies by roughly a factor of 3 from the top of the pedestal to the bottom, with the profile of 
dimensionless n = ni/n\{r = —00) shown in[8la. This density profile implies the electric field 
profile shown in figure Ob, which reaches a minimum of ~ —0.5viBg/c in the pedestal center. 
The coUisionality ranges from 0.5 — 0.15 over the domain. (We choose this arbitrary range 
close to one just to emphasize that is not formally large or small in this formulation.) In 
these plots, the radial coordinate r/pg is defined by r/pg = ZeBQ{miCViI)~^ip where Bq is 
the toroidal field on axis. The radial location r = is an arbitrary minor radius, (here the 
middle of the pedestal), not the magnetic axis. The sink used is fl5^ with 7 = O.lut, where 
Ut = fiV||^ is the transit frequency. The simulation is nearly converged by t = 30/ciJt, but 
very small changes in the results continue until t = 200/(jJt. We plot results for t = 200/a;t 
since doubling this duration produces no visible change to the results. By t = 200/a;t, the 
residual, which we define as a sum over all phase-space grid points of \dg/dt\, has been 
reduced to 0.05% of its initial value. 
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It is also important to verify that the code has converged with respect to the many other 
numerical parameters. Figure Eld shows the parallel flow coefficient k\\ at the outboard mid- 
plane for 11 global runs, all with the same physics parameters, but varying each numerical 
parameter by a factor of two: simulation duration (tmax)? time step (dt), artificial radial vis- 
cosity, number of poloidal modes {N6), number of Legendre polynomials in the Rosenbluth 
potentials (NL), number of grid points in v, ^, and ip {Nv, NC,, and Nip), and domain size 
in speed (fmax) and radius (V^max)- The changes are barely perceptible, demonstrating very 
good convergence. For comparison, the profile computed by the local code is also plotted, 
calculated by solving (i.e. a single linear system solve) at each radial grid point. The 
local coefficient varies across the pedestal due to the change in collisionality. Resolution pa- 
rameters were, unless doubled, N9 =6, Nip =29, A^,^ =25, Nv =16, r^ax = 4:p9, v^^^ = 5vi, 
NL =2, and dt = O.Olwt. Running in Matlab on a single Dell Precision laptop with In- 
tel Core 17-2860 2.50 GHz CPU and 16 GB memory, the base case global simulation took 
roughly 3 hours to reach t = 200aJt, though runs could undoubtedly be greatly expedited 
if the code were parallelized and rewritten in Fortran. Work to this end is underway. The 
local solver for these parameters took 0.5 seconds per radial grid point. 

Figures [Hld-f show the heat fiux kgn"^ and the fiow coefficients k\\ and kg at the outboard 
and inboard midplanes. (It is radial variation in the heat flux kg-n^ and not kg itself that 
determines local heating, as shown by ([5]) and (l53l)). Outside of the pedestal, as expected, 
these coefficients agree with the local prediction, and k\\ and kg are equal and poloidally 
invariant. In the pedestal, however, all coefficients are substantially altered from the local 
prediction; k\\ and kg differ and vary poloidally. The radial heat flux proflle is flattened 
relative to the local prediction. 

For the parameters used here, k\\ and kg change sign in the simulation within the pedestal. 
These coefficients may have either sign in conventional theory, as shown in figure [21 depend- 
ing on collisionality and magnetic geometry. In both the local and global cases, the integrand 
v\\g that determines k\\ oc J d^vv\\g is positive in part of phase space and negative elsewhere, 
and the balance between these regions determines the overall sign. 

Structure with a radial scale comparable to pg is observed in the flow coefficients. Ions 
"communicate" over distances comparable to the orbit width ~ y/^pg, and so the effects of 
the driving electric field well are felt outside of the well itself, with influence decaying on the 
orbit width scale. The behavior of the flow coefficients on either side of the well need not 



29 





1 


a 


0.5 








-0.5 




0.5 


> 






n 
u 




0.2 


II 




II 

as 


U 










-0.2 




1 




8 
U.o 




0.6 


c 


a" 

M 


U.4 




0.2 









0.4 


CD 


0.2 


O 







-0.2 




-0.4 




-4 




a) 



b) 



c) 




f) 



■ Local result 

■ Global solution 
2xNe 
2x N\|/ 

■ 2xN^ 

■ 2xNv 

■ 2x NV|/ and ymax 
2x Nv and vniax 
2x viscosity 
2x NL 
2x tmax 
Half dt 





— Local k = k 


Global k at 6=0 

9 




II s 

— Global \ at 6=0 


Global at 6=71 




— Global k|| at 6=7r 



-2 







FIG. 8. (Color online) a) Equilibrium density profile for the global calculation, normalized to 
its value at the left boundary, b) Normalized radial electric field —cI{v\B{^)^^ d^^^j d"^ . c) Profile 
of Vtf. d) Parallel flow coefficient evaluated at the outboard midplane {0 = 0). Black dashed 
curve indicates the result of the local neoclassical code. The other curves (nearly indistinguishable) 
demonstrate the convergence of the global code to the various numerical resolution parameters, e) 
Radial heat flux, with the same legend as d). f) Normalized poloidal flow ke and evaluated at 
the outboard and inboard {9 = vr) midplanes. 



be monotonic, as the mean flow adjacent to the well arises from a complicated interplay of 
particles entering from regions of differing coUisionality, some particles directly affected by 
the electric field and some not. 

The flow coefficients are also observed to be non-monotonic functions of r. This behavior 
is not unreasonable given the radial localization of the electric field well : even the local 
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FIG. 9. (Color online) Poloidal variation of the parallel flow coefficient and normalized poloidal 
flow ko at two radial locations straddling the pedestal. 

distribution function is a non-monotonic function of r, since it is a complicated function 
of the radially varying coUisionality, so in the global case the flow coefficients need not be 
monotonic in r. 

Figure [9] shows the poloidal variation of the flow coefficients near the pedestal top and 
bottom. As discussed in the appendix, the drift terms in the kinetic equation break the sym- 
metry which the distribution function possesses in the local case, and so the flow coefficients 
need not be even or odd in 6. 

The various components of the mass conservation equation (1261) were each independently 
computed from g: V ■ / d?vv\\gb, V ■ / d^WEog, V ■ / d^vv^g, and / d^vS. The first 
three of these integrals summed to nearly zero everywhere in space, with the sink integral 
negligible in magnitude compared to the others. Thus, the sink has little effect on the mass 
conservation relation that effectively determines the flows. The f m integral was intermediate 
in magnitude, leaving a dominant balance between the V- / d^vv\\gb and V- / d^vv^og terms. 
The density J d^v g has a oc cos(6') behavior in the pedestal, resulting in V ■ / d^vvEog oc 
[d/dO) J d^v g oc — sin(6'). To balance this term in the mass conservation law, k\\ must 
develop a — cos(6') structure, which can be seen in flgure ([9]). Outside of the pedestal, the 
veo term becomes negligible due to the reduced \d(^Q/dip\, so this drive for poloidal variation 
in k\\ is absent. 

Figure [TD] shows how the results are altered when different choices are made for the sink. 
For the sink we show results for 7 = O.lwt (the value used for all other plots) and for 
7 = cjf We also show results for the alternative sink ( 155|) . For comparison, results are also 
shown for a run in which no sink was included. For this run the code did not converge in 
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FIG. 10. (Color online) Using the same equilibrium density profile (a) and Er profile (b) as in 
Figure ([8]), the parallel flow coefficients (c-d) and heat flux (e) show some minor dependence on 
the choice of sink. However, qualitative features such as the well in k^^ at the outboard midplane 
are robust. 

time, due to the constant loss of heat described in section IVIIIt so it was stopped at t = SOcut 
(a time before the heat loss becomes excessive, but after the runs with sinks have nearly 
converged.) The various options yield results that show the same qualitative modification 
of the coefficients: a well develops in k\\ and kg at the outboard midplane, and the heat flux 
profile is fiattened relative to the local prediction. 

As a further test of the code, we repeated the numerical calculation, treating /i as the 
unknown quantity instead of (? , in which case the inhomogeneous term in the equation is 
—V]j- V/Mi instead of Ci(F). Despite the very different phase-space structure of these two 
source terms, the numerical results from the two approaches agreed, as they should. 
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X. DISCUSSION 



In this work we have demonstrated a method to extend neoclassical calculations to in- 
corporate finite-orbit-width effects in a transport barrier for the case ol Bq <^ B (which is a 
good approximation in standard tokamaks) and a relatively weak ion temperature gradient. 
The method is implemented in a new continuum 5f code. Operator splitting is used to im- 
prove numerical efficiency, and we have demonstrated that excellent convergence is feasible 
for experimentally relevant parameters. By construction, the method exactly reproduces 
conventional (local) results in the appropriate limit of weak radial gradients. The Rosen- 
bluth potentials are solved for along with the distribution function at each step, allowing 
use of the full linearized Fokker-Planck-Landau collision operator. 

A principal finding of this work is that the parallel and poloidal fiows may differ signif- 
icantly from the conventional predictions. While the coefficients of the poloidal fiow and 
(iTi/rf-j/^-driven parallel fiow are equal in conventional theory, in the pedestal these two coef- 
ficients {k\\ and kg) differ. And, while the poloidal variation of the poloidal ffow is oc Bg in 
the core, the same is not necessarily true in the pedestal. The poloidal variation of the ffow 
is effectively determined by mass conservation, and in the pedestal, two new terms become 
important which are normally neglected: E x B convection of the poloidally varying den- 
sity, and radial variation of the particle ffux (which can be related to diamagnetic ffow from 
the correction to the pressure). These effects cause the parallel ffow coefficient k\\ to take 
a well-shaped radial profile, with different magnitude and (for some parameters) opposite 
sign, relative to the conventional local neoclassical result. In addition, the ffow coefficients 
exhibit a strong poloidal variation not previously found. While this poloidal variation re- 
sembles cos(6'), it contains other harmonics and is asymmetric about the midplane due to 
the magnetic drift. 

These issues may be important for comparisons of experimental pedestal fiows to 
theory^i^. In general, the fiow coefficients may differ in both magnitude and sign rel- 
ative to local theory, as shown in figure [HI The fiuid fiow exhibits strong shear, with radial 
variation on the pg scale. 

Associated with the modification to the fiow, the parallel current is also modified. Due to 
the additional terms which must be included in the mass conservation equation, the usual 
division of the parallel current into Pfirsch-Schliiter and Ohmic-bootstrap components is 
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modified, as shown in Eq. (jH]). In addition, the dTi/dip contribution to the bootstrap 
current is altered, as shown in P5|) . In the 5f formulation, the ion temperature scale length 
cannot be as small as the density scale length in the pedestal, so these modifications to the 
parallel current are modest. However, similar changes to the current would presumably occur 
in a full-/ calculation when tt ~ pe-, giving order-unity changes to the Pfirsch-Schliiter and 
bootstrap currents in that case. This issue needs to be examined in future studies. 

In the development of this work, the local code was also used to test several analytic 
expressions for conventional neoclassical theory. The flow and heat flux coefficients k\\ = 
kg = l.n,kq = 0.66 derived using the momentum-conserving pitch-angle scattering model 
for collisions are a poor approximation unless e is ^ 0.1. Expressions f|T6|) -f lT7j) are a much 
better approximation at realistic aspect ratio. The Chang-Hinton heat flux captures the 
trends at finite e and i/^, well and gives results correct to within 20%, at least for the circular 
concentric flux surface model. The semi-analytic local formulae of Sauter et al^ for the flow 
and bootstrap current coefficients were found to be in excellent agreement with our local 
code when < 0.3, but some disagreement was found for 0.3 < i^* < 100. For pedestal 
profiles typical of DIII-D, the Sauter bootstrap current formula closely agreed with our code 
at low collisionality, i^*e < 1? but the Sauter formula can give a bootstrap current more than 
twice ours when > 10. The Sauter formulae are intended to reproduce results from a 
code based on the same physical model as our conventional local code. 

There are many ways in which the global calculations can be extended. First, it would be 
useful to include impurities, for it is typically the impurity flow that is measured rather than 
that of the main ions, and the the flows of different ion species may be significantly different^. 
Also, the presence of impurities can introduce a direct density gradient dependence^ to g. 
Second, the method should be extended to allow strong temperature gradients {r^ ~ pe)- 
Doing so will require the full bilinear collision operator and a full-/ treatment. However, as 
the weak- 7]' case is less complicated to analyze, due to the linearity of the kinetic equation, 
thorough understanding of this limit using the present approach is important for benchmark- 
ing future more sophisticated full-/ codes. Finally, studies of the velocity-space structure 
responsible for fluxes in turbulence codes may yield more accurate forms of the sink term 
needed in our approach. 
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Appendix A: Null space and symmetry of the distribution 

The local drift-kinetic equations ([1]) and ([3]) have two null solutions /mi and v'^fui, mean- 
ing that the discretized matrix for the local code should be nearly singular. In practice the 
matrix is still sufficiently well conditioned that the linear system may be solved without 
a problem, yielding a distribution function that contains a small amount of the two null 
solutions. For many applications this may not be a concern, because these null solutions do 
not contribute to the heat flux and flow. 

For an up-down symmetric tokamak (i.e. if B and \/\\6 are both unchanged under 6 — )■ 
—0), then a symmetry exists in the local kinetic equations: if fi{6,C,,v) is a solution, then 
so is —fi{—6, v) (and similarly for g). This property can be exploited to simultaneously 
eliminate the null space from the matrix and to reduce its size^. This is done by forcing /i 
to have the above symmetry by representing it as a sum of two types of modes: those that 
are even in 6 and odd in ^, and those that are odd in 6 and even in ^. The two null solutions 
do not possess this symmetry, so they are automatically excluded. Furthermore, the matrix 
size is reduced without loss of resolution. For example, the C, grid can be reduced to only 
cover the interval [0, 1] instead of [—1, 1] if all sin(M^) and cos(M^) modes are retained. 
The odd-^ (i.e. sin(M^)) modes are forced to be even in ^ by application of the boundary 
condition dfi/dC, = at ,^ = 0, and the even-6' (i.e. cos(M6')) modes are forced to be odd in 
^ by application of the boundary condition /i = at ^ = 0. 

Even if the parity of the solution is not enforced automatically by the discretization in 
this manner, the null solutions can still be excluded by enforcing parity as follows. Given 
a numerical solution fi{6,$,,v) that contains some of the null solutions, the combination 
[fi{9, ^, v) — fi{—0, v)] /2 can be formed; the result will also satisfy the kinetic equation 
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but have the desired parity. 

In the global case, the symmetry of the kinetic equation is broken by the drift terms. 
However, the local operator has no null space in the initial-value-problem formulation due 
to the extra contribution on the matrix diagonal from the time derivative. 



Appendix B: Conservation laws for the global drift-kinetic equation 

Here we sketch the derivation of general conservation equations, from which fl2B]) . (13^ . 
and fl53|l can be obtained. For most of this appendix we do not assume axisymmetry, we 
retain radial variation of magnetic quantities, and we do not require B ■ V$ = 0. We 
do require the electric field to be electrostatic and we assume d^/dt and dB/dt can be 
neglected. The derivation applies both to the full-/ and 6f contexts, since the necessary 
integrals of both the bilinear and linearized ion-ion collision operators vanish. 

We begin with the ion drift-kinetic equation 

df,/dt + V||/i + «d ■ V/i = + S (Bl) 

where = {v\\/Q)'V x {v\\b) and C is either the bilinear or linearized Fokker- Planck- 
Landau operator. Gradients are all performed at fixed fi and total energy W = mit;^/2 + Ze$ 
(including the total potential $, not just $o)) so i>d includes both the magnetic drift 
and E X B drift ve- This form of is convenient because it makes the kinetic equation 
conservative without cumbersome higher-order terms. This includes an incorrect 0(/3) <C 
1 parallel magnetic drift, but this component of the magnetic drift is typically unimportant 
compared to parallel streaming motion, and in fact ■ V6 is precisely zero in the model 
magnetic geometry we use in the code. It is convenient to first rewrite 

^||V||/i + V, . V/i = |V ■ [f,B + X V/i) . (B2) 

Then J d^v is applied to f IBip . annihilating C. Notice 

27r ^ p{W-Ze<f)/B o 

^ JZe'}> Jo ^11 

where a = sgn(f||). The divergence in f lB2p may be pulled in front of the integrals in (lB3p . 
as the contributions from differentiating the integration limits all vanish either due to the a 
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sum or because f|| = at the lower limit of W. Application of several vector identities to 
the B X V/i term then yields a mass conservation equation: 

^ (/ '^'''^0 (/ = / (^4) 

Flux surface averaging and neglect of $i then gives ( l52l) . 

An energy conservation equation may be obtained by observing that the above derivation 
of the mass conservation equation is essentially unchanged if J d?v W is applied to f IBip in 
place of / dPv. Subtracting Ze'l'x ( lB4l) from the result, one obtains 

!(/ ^'^^^^)+^{/ dMv\\h + v,]^f^ + (^j dMv\\h + v,]f^.ZeV^ = I dh^S. 

(B5) 

In the special case of axisymmetry and V||$ = 0, flux surface averaging and integration in 
ip then gives ( l53l) . 

To obtain the momentum conservation equation, it is convenient to specialize to ax- 
isymmetry at the start, taking the J d^v{Iv\\/ B) moment of (IBip . and using ■ V/i = 
(f||/i?)V ■ [/i(mic/Ze)V x (f||b)] instead of (IB2|) . The divergence may be brought in front 
of the W and /i integrals as before. Noting ■ 'S/{Iv\\/ B) = and v\\\/\\{Iv\\/Q) = ■ V^/j, 
the result may be written 

I (/ ^'^^^^)+^-(/ + / ^'^/i«d-V^ = I dh^-^S. (B6) 

This result holds in axisymmetry even if V||$ and/or V||/ are nonzero. 



Appendix C: Convenient gauge 



Here we prove that the gauge may always be chosen so 



^ll = ^(^l|5>-V||$. (CI) 



Axisymmetry is not required, and the loop voltage need not be uniform. The utility of ( !CT 
is that the inductive part of E\\ has simple spatial variation oc B. 
Suppose we begin in a different gauge, denoted by tildes, in which 



E = ~c-^dA/dt-Vi\^. 
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(C2) 



We may transform to a new gauge using $ = ^ — dx/dt and A = A — Vx ioT a generator 
X- We choose 



where the integrand is evaluated at t' and 6' rather than t and 6. We must verify fICSp is 
single- valued in 6* so $ is single- valued. To this end, notice {B ■ ( )) applied to f lC2|) gives 
(^E\\B'^ = — (^B ■ dA/dtJ. Therefore x{^ = 27r) = = x{^)^ so x is indeed periodic. 
Applying V|| to $ = <l - 9x/9t with S-(IC2]) and then gives ([Cl]) as desired. 
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